%% Producing Figure 5: Impulse response to show procyclical leverage

close all
warning off

load('no_policy')

%% Set up variables

options = optimset;

before = 4;
after  = 20;

shocks_run           = zeros(1,2+before+after);
shocks_run(2+before) = 1;

kv   = zeros(1,2+before+after);
rbv  = zeros(1,2+before+after);
nv   = zeros(1,2+before+after);
xiv  = zeros(1,2+before+after);
hv   = zeros(1,2+before+after);
cv   = zeros(1,2+before+after);
yv   = zeros(1,2+before+after);
iv   = zeros(1,2+before+after);
lv   = zeros(1,2+before+after);
pibv = zeros(1,2+before+after);
runv = zeros(1,2+before+after);
rhseev = zeros(1,2+before+after);
rkhat  = zeros(1,2+before+after);
rkhats = zeros(1,2+before+after);
prv  = zeros(1,2+before+after);
av   = zeros(1,2+before+after);
erk  = zeros(1,2+before+after);
z_temp = zeros(1,2+before+after);
spreadv = zeros(1,2+before+after);

x    = zeros(2+before+after,npol);
x_nonlog    = zeros(2+before+after,npol);
x_allnonlog = zeros(2+before+after,npol);

%% Simulation

kv(1)  = stss_k;
rbv(1) = stss_rb;
nv(1)  = stss_n;

lv(1)  = kv(1)/nv(1);
av(1)  = 0;

    for t=2:2+before+after
        av(t) = rho_a*av(t-1) + sigma_a*shocks_run(t);
        hv(t)  = ((1-alpha)/psi*exp(av(t))*kv(t-1)^(alpha))^(1/(alpha+1/nu));
        yv(t)  = exp(av(t)) * kv(t-1)^(alpha)*hv(t)^(1-alpha);

        rkhats(t) = log ( rbv(t-1) * (1-1/lv(t-1)) * ( 1+lambda*(1-gamma) ) - (1-delta) );
        rkhat(t)  = log ( alpha*((1-alpha)/psi)^((1-alpha)/(alpha+1/nu)) ) ...
                  + (1+nu)/(alpha*nu+1)*av(t) - (1-alpha)/(alpha*nu+1)*log(kv(t-1));
        xiv(t)  = exp(rkhat(t))/exp(rkhats(t));

        K_nonlog  = [log(kv(t-1)) rbv(t-1) log(nv(t-1)) xiv(t)];
        x0_nonlog = (K_nonlog).^expmx;
        x_nonlog(t,:)  = prod(x0_nonlog,2)';
        
        K_allnonlog  = [kv(t-1) rbv(t-1) nv(t-1) xiv(t)];
        x0_allnonlog = (K_allnonlog).^expmx;
        x_allnonlog(t,:)  = prod(x0_allnonlog,2)';

        K  = [kv(t-1) rbv(t-1) nv(t-1) xiv(t)];
        x0 = (log(K)).^expmx;
        x(t,:)  = prod(x0,2)';
        if xiv(t) >= 1
            pibv(t)  = ( exp(rkhat(t)) + 1-delta )*lv(t-1)*nv(t-1) - rbv(t-1)*(lv(t-1)-1)*nv(t-1) - nv(t-1);

            if pibv(t) > 0
                rhseev(t) = exp(x(t,:)*solved_eta1_rhsee);
                rbv(t)    = exp(x(t,:)*solved_eta1_rbar);
                nv(t) = chi1*( chi0*pibv(t) + nv(t-1) ) + n0;
            else
                rhseev(t) = exp(x(t,:)*solved_eta3_rhsee);
                rbv(t)    = exp(x(t,:)*solved_eta3_rbar);
                nv(t) = chi1*(      pibv(t) + nv(t-1) ) + n0;
            end

        else
            pibv(t)  = ( exp(rkhat(t)) + 1-delta )*lv(t-1)*nv(t-1) - rbv(t-1)*(1+lambda*(1-gamma))*(lv(t-1)-1)*nv(t-1) - nv(t-1);
            nv(t) = n_ss*(1 - nd/100);
            rhseev(t) = exp(x(t,:)*solved_eta2_rhsee);
            rbv(t)    = exp(x(t,:)*solved_eta2_rbar);
            runv(t) = 1;
        end

        cv(t) = rhseev(t)^(-1) + psi*hv(t)^(1+1/nu)/(1+1/nu);
        iv(t) = yv(t) - cv(t);
        kv(t) = (1-delta)*kv(t-1) + iv(t);

        lv(t)  = kv(t)/nv(t);

        erk(t) = log ( alpha*((1-alpha)/psi)^((1-alpha)/(alpha+1/nu)) ) ...
               + (1+nu)/(alpha*nu+1)*rho_a*av(t) - (1-alpha)/(alpha*nu+1)*log(kv(t));

        rbv(t) = fzero('func_Rb',rbv(t),options,lv(t),erk(t),lambda,gamma,delta,sigma_rk);   % deposit interest rate

        rkhats_next = log ( rbv(t) * (1-1/lv(t)) * ( 1+lambda*(1-gamma) ) - (1-delta) );
        z_temp(t)     = (rkhats_next - erk(t)) / sigma_rk;
        cdf = normcdf(z_temp(t));
        prv(t) = cdf;
        spreadv(t) = exp(erk(t)) - rbv(t);

    end
    % End of T period simulation

%% plot resulsts

close all

figure('Name','investment','Units','Inches','OuterPosition',[7 0 8 8])
plot(-before:after,iv(2:end)/stss_i*100-100,'Linewidth',4)
hold on
plot(-before:after,zeros(1,1+before+after),'k:','Linewidth',2);
hold on
plot([0 0],[-1 4],'k-.','Linewidth',2);
set(gca,'FontSize',40);
xlim([-before after])
xticks(-before:4:after)
xtickangle(0)
ylim([-1 4]);
yticks(-1:1:4)
grid on
ytickformat('percentage')

figure('Name','net worth','Units','Inches','OuterPosition',[7 0 8 8])
plot(-before:after,nv(2:end)/stss_n*100-100,'Linewidth',4)
hold on
plot(-before:after,zeros(1,1+before+after),'k:','Linewidth',2);
hold on
plot([0 0],[-0.2 0.8],'k-.','Linewidth',2);
set(gca,'FontSize',40);
xlim([-before after])
xticks(-before:4:after)
xtickangle(0)
ylim([-0.2 0.8]);
yticks(-0.2:0.2:0.8)
grid on
ytickformat('percentage')

figure('Name','leverage','Units','Inches','OuterPosition',[7 0 8 8])
plot(-before:after,lv(2:end),'Linewidth',4)
hold on
plot([0 0],[10 10.4],'k-.','Linewidth',2);
set(gca,'FontSize',40);
xlim([-before after])
xticks(-before:4:after)
xtickangle(0)
ylim([10 10.4]);
yticks(10:0.1:10.4)
grid on

figure('Name','bank return','Units','Inches','OuterPosition',[7 0 8 8])
plot(-before:after,rbv(2:end),'Linewidth',4)
hold on
plot([0 0],[1.005 1.025],'k-.','Linewidth',2);
set(gca,'FontSize',40);
xlim([-before after])
xticks(-before:4:after)
xtickangle(0)
ylim([1.005 1.025]);
yticks(1.005:0.005:1.025)
grid on

figure('Name','run probability','Units','Inches','OuterPosition',[7 0 8 8])
plot(-before:after,prv(2:end)*100,'Linewidth',4)
hold on
plot([0 0],[1.2 1.6],'k-.','Linewidth',2);
set(gca,'FontSize',40);
xlim([-before after])
xticks(-before:4:after)
xtickangle(0)
ylim([1.2 1.6]);
yticks(1.2:0.1:1.6)
grid on
ytickformat('percentage')
%title('output','Fontsize',40)

figure('Name','TFP','Units','Inches','OuterPosition',[7 0 8 8])
plot(-before:after,av(2:end)*100,'Linewidth',4)
hold on
plot(-before:after,zeros(1,1+before+after),'k:','Linewidth',2);
hold on
plot([0 0],[-0.2 1.2],'k-.','Linewidth',2);
set(gca,'FontSize',40);
xlim([-before after])
xticks(-before:4:after)
xtickangle(0)
ylim([-0.2 1.2]);
yticks(-0.2:0.2:1.2)
grid on
ytickformat('percentage')
